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We apply the "weighted ensemble" (WE) simulation strategy, previously employed in the context of molecular 
dynamics simulations, to a series of systems-biology models that range in complexity from one-dimensional to 
a system with 354 species and 3680 reactions. WE is relatively easy to implement, does not require extensive 
hand-tuning of parameters, does not depend on the details of the simulation algorithm, and can facilitate the 
simulation of extremely rare events. 

For the coupled stochastic reaction systems we study, WE is able to produce accurate and efficient approxi- 
mations of the joint probability distribution for all chemical species for all time t. WE is also able to efficiently 
extract mean first passage times for the systems, via the construction of a steady-state condition with feedback. 
In all cases studied here, WE results agree with independent calculations, but significantly enhance the precision 
with which rare or slow processes can be characterized. Speedups over "brute-force" in sampling rare events via 
the Gillespie direct Stochastic Simulation Algorithm range from ~10 12 to ~10 20 for rare states in a distribution, 
and ~10 2 to ~10 4 for finding mean first passage times. 



I. Introduction 

Stochastic behavior is an essential facet of biological pro- 
cesses such as gene expression, protein expression, and epi- 
genetic processes IfTl TBl . Stochastic chemical kinetics sim- 
ulations are often used to study systems biology models of 
such processes fTsHTTIl . One of the more common stochas- 
tic approaches, and the one employed in the present study, is 
the stochastic simulation algorithm (SSA), also known as the 
Gillespie algorithm lfl5l [181 [191 . 

As stochastic systems biology models approach the true 
complexity of the systems being modeled, it quickly be- 
comes intractable to investigate rare behaviors using naive 
("brute-force") simulation approaches. By their very nature, 
rare events occur infrequently; confoundingly, rare events 
are often those of most interest. For example, the switch- 
ing of a bistable system from one state to another may hap- 
pen so infrequently that running a stochastic simulation long 
enough to see transitions is (extremely) computationally pro- 
hibitive [20|. This impediment only grows as model complex- 
ity increases, and as such it poses a serious hurdle for systems 
models as they grow more intricate. 

Several approaches to speeding up the simulation of rare 
events in stochastic chemical kinetic systems exist. A vari- 
ety of "leaping" methods can, by taking advantage of approx- 
imate time-scale separation, accelerate the SSA itself ET1 - 
[28). Kuwahara and Mura's weighted stochastic simulation 
(wSSA) method [29 1 was refined by Gillespie and Petzold et 
al. E0-33], and is based on importance sampling. The for- 
ward flux sampling method of ten Wolde et al. Il20l [3414361 
uses a series of interfaces in state-space to reduce computa- 
tional effort, as does the non-equilibrium umbrella sampling 
approach 1371 l38l . 

Rare event sampling is also an active topic in the field of 
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molecular dynamics simulations, and many approaches have 
been proposed. Of the approaches that do not irreversibly 
modify the free energy landscape of the system, some no- 
table methods include dynamic importance sampling [39], 
milestoning Pol , transition path sampling 1411 . transition 
interface sampling [42J, forward flux sampling |36|, non- 
equilibrium umbrella sampling ll38ll . and weighted ensemble 
sampling [ 43-50 1. For a summary of these methods, see BTl . 
Many of the ideas behind these techniques are not exclusive to 
molecular dynamics simulations, and can be adapted to study- 
ing stochastic chemical kinetic models. For example, dynamic 
importance sampling seems to be closely related to wSSA. 

Because of its relative simplicity and potential simplicity 
in sampling rare events, we apply one of these methods, the 
weighted ensemble algorithm (WE) to well-established model 
systems of stochastic kinetic chemical reactions. These mod- 
els range in complexity from one species and two reactions, to 
354 species and 3680 reactions. For the systems studied, WE 
proves many orders of magnitude faster than SSA simulation 
alone, offers linear parallel scaling, returns full distributions 
of desired species at arbitrary times, and can yield mean first 
passage times (MFPTs) via the setup of a feedback steady- 
state. 



II. Methodology 

The methods employed are described immediately below, 
while the models are specified in Sec. [Hi] 

A. Stochastic Chemical Kinetics & BioNetGen 

Stochastic chemical kinetics occupies a middle-ground in 
the realm of chemical simulation, between very explicit, and 
costly, molecular dynamics (MD) simulations and the deter- 
ministic formalism of reaction rate equations (RRE). Stochas- 
tic chemical kinetics attempts to account for the randomness 
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inherent in chemical reactions, without trying to explicitly 
model the spatial structure of the reacting species. It is many 
orders of magnitude faster than MD simulations, but much 
slower than the RRE approach. It is an ideal method to use for 
modeling the effects of low concentrations (or copy numbers) 
of chemical reactants, while ignoring the effects of specific 
spatial distribution. 

Stochastic chemical kinetics models can be solved exactly 
for sufficiently simple systems using the Chemical Master 
Equation (CME), and approximately (for all systems) using 
Gillespie's direct stochastic simulation algorithm (SSA) lfl5l 
COD [19] . The SSA samples the CME exact solution by mod- 
eling stochastic chemical kinetics in a straightforward man- 
ner, and yields trajectories of species concentrations that con- 
verge to the RRE method in the limit of large amounts of re- 
actants. In brief, the SSA iteratively and stochastically deter- 
mines which reaction fires when reactions occur, by sampling 
from the exponential distribution of waiting times between re- 
actions. For a detailed explanation of the SSA, see ifTSl . 

We employ the rule-based modeling and simulation pack- 
age BioNetGen |52| to simulate both our toy and complex 
models. Rule-based modeling languages allow the specifi- 
cation of biochemical networks based on molecular interac- 
tions. Rules that describe those interactions can be used to 
generate a reaction network that can be simulated either as 
RREs or using the SSA, or the rules can be used directly to 
drive stochastic chemical kinetics simulations. BioNetGen 
has been applied to a variety of systems, such as the aggre- 
gation of membrane proteins by cytosolic cross-linkers in the 
LAT-Grb2-SOS 1 system l53l . the single-cell quantification of 
IL-2 response by effector and regulatory T cells 11541 . the anal- 
ysis and verification of the HMGB1 signaling pathway Il55l , 
the role of scaffold number in yeast signaling systems ll56l . 
and the analysis of the roles of Lyn and Fyn in early events 
in B cell antigen receptor signaling [57|. We employ BioNet- 
Gen's implementation of the direct SSA to propagate the dy- 
namics in our systems. 



B. Weighted Ensemble (WE) 

WE is a general-purpose protocol used in molecular dy- 
namics simulations [ 44-46 48 - 50] that we adapt here to the 
efficient sampling of dynamics generated by chemical kinetic 
models. In brief, WE employs a strategy of "statistical natural 
selection" using quasi-independent parallel simulations which 
are coupled by the intermittent exchange of information. The 
intermittency leads directly to linear parallel scaling. Impor- 
tantly, the simulations are coupled via configuration space (es- 
sentially the "phase space" of the system in physics language 
or the "state-space" in cell and population modeling). This 
type of coupling permits both efficiency and a large degree 
of scale independence. The efficiency results from distribut- 
ing trajectories to typically under-sampled parts of the space, 
while scale independence is afforded because every type of 
system has a configuration or state-space. 

WE's strategy of statistical natural selection or statistical 
ratcheting is schematized in Fig. [T] First, the space is divid- 



ed/classified into non-overlapping "bins" which are typically 
static, although dynamic and adaptive tessellations are possi- 
ble 1351 . A target number of trajectories, M talg , is set for each 
bin. Multiple trajectories are initiated and each is assigned a 
weight so that the sum of weights is one. Trajectories are then 
simulated independently according to the desired dynamics 
(e.g., molecular dynamics or SSA) and checked intermittently 
(every t units of time) for their location. If a trajectory of 
weight w is found to occupy a previously unoccupied bin, that 
trajectory is replicated to obtain the target number of copies, 
M targ , for the bin. Daughter trajectories' weights are set to 
w/M talg , to sum to the weight of the parent trajectory. If a 
bin is occupied by more than the target number, trajectories 
must be pruned in a statistical fashion maintaining the sum 
of weights. Specifically, the two lowest weight trajectories 
are "merged" by randomly selecting one of them to survive, 
with probability proportional to their weights, and the surviv- 
ing trajectory absorbs the weight of the pruned one. This pro- 
cess is repeated as needed, and maintains an exact statistical 
representation of the evolving distribution of trajectories 11451 . 

Setting up a WE simulation requires selection of state- 
space binning, trajectory multiplicity, and timing parameters. 
In our simulations, we chose to divide the state-space of an 
iV-dimensional system into one- or two-dimensional regular 
grids of non-overlapping bins. It is possible to use non- 
Cartesian bins, and to adaptively change the bins during sim- 
ulation [45 48 1, but for simplicity we did not pursue any such 
optimization. Specific parameter choices for each model are 
given in Sec.|Hl| 

The weighted ensemble algorithm can be outlined fairly 
concisely. Let M targ be the target number of segments in each 
bin, A'bins the number of bins, whose geometry are defined by 
the grid G gl -id, r the time-step of an iteration of WE, and Alters 
the total number of iterations of WE. The WE procedure also 
requires an initial state of the system, Xq, which in our case is 
a list of the concentrations of all the chemical species in the 
system. 

procedure WE(Mters, t, G grid , M targ , x ) 
for i = 1 . . . Af itels do 

for each populated bin in G gl -id do 

propagate dynamics for all trajectories 

for each bin in G gl -id do 

if bin population = or M talg then 

do nothing 
else if bin population < M tare then 

replicate trajectories until bin pop. = M talg 
maintain sum of weights in each bin 
else if bin population > M targ then 

merge trajectories until bin pop. = M targ 
maintain sum of weights in each bin 
save coordinates and weights of each trajectory 
return trajectory coordinates & weights for each iter. 

The replicating and merging of trajectories in the above al- 
gorithm are done randomly, according to the weight of each 
trajectory segment in a given bin, which has been shown not 
to bias the dynamics of the ensemble ll43l l46l . 

When WE is used to manage an ensemble of trajectories, 
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FIG. 1. Weighted ensemble (WE) simulation depicted for a configuration/state-space divided into bins. Multiple trajectories are run using any 
dynamics software (here we use the SSA in BioNetGen) and checked every r for bin location. Trajectories are assigned weights (symbols - 
see legend) that sum to one and are split and combined according to statistical rules that preserve unbiased kinetic behavior. 



there are two time-scales of immediate concern: the period 
at which trajectory coordinates are saved, and the period r at 
which ensemble operations are performed. These two time- 
scales can be different, but for simplicity we set them to be 
the same, and select t such that it is greater than the average 
event firing rate for the SSA. When we refer to the time-step, 
or iteration of a process, we are referring to the t of Fig.[T] 

WE can be employed in a variety of modes to address dif- 
ferent questions. Originally developed to monitor the time 
evolution of arbitrary initial probability distributions |43l , i.e. 
non-stationary non-equilibrium systems, WE was generalized 
to efficiently simulate both equilibrium and non-equilibrium 
steady-states |46l . In steady-state mode, mean first passage 
times (MFPTs) can be estimated rapidly based on simula- 
tions much shorter than the MFPT using a simple rigorous 
relation between the flux and MFPT |46l . Steady-states can 
be attained rapidly, avoiding long relaxation times, by using 
the inter-bin rates computed during a simulation to estimate 
bin probabilities appropriate to the desired steady-state; tra- 
jectories are then reweighted to conform to the steady-state 
bin probabilities |46|. Both of these methods are described in 
more detail below. 



1. Basic WE: Probability Distribution Evolving in Time 

Perhaps the simplest use of a weighted ensemble of tra- 
jectories is to better sample rare states as a system evolves 
in time, specifically the states corresponding to extreme val- 
ues of the binning coordinate. The SSA itself samples the 



exact distribution, but its sampling is concentrated about the 
mode(s) of the distribution. The SSA naturally - and correctly 
- samples rare states infrequently. By using WE to split up the 
state-space, however, one can resample the distribution at ev- 
ery time step t, selecting those trajectories that advance along 
a progress coordinate for more detailed study, but doing so 
without applying any forces or biasing the trajectories or the 
distribution. Essentially, WE appropriates much of the effort 
that brute-force SSA devotes to sampling the central compo- 
nent of the distribution, repurposing it to obtain better esti- 
mates of the tails. 

This basic use of WE requires none of the "tricks" we apply 
in later sections, such as using reweighting techniques to ac- 
celerate obtaining a steady-state. We apply basic WE to some 
of our systems - particularly, but not exclusively, to those that 
are not bistable. 



2. Steady-State 

The mean first passage time (MFPT) from state A to state 
B is a key observable. It is equal to the inverse of the flux 
(of probability density) from state A to state B in steady-state 

El, 



MFPT A ^ B 
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(1) 



This relation provides the weighted ensemble approach the 
ability to calculate MFPTs in a straightforward manner. Dur- 
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ing a WE run, when any trajectories (and their associated 
weights) reach a designated target area of state-space (or 
"state fi"), they are removed and placed back in the initial 
state ("state A"). Eventually, such a process will result in a 
steady-state flow of probability from state A to state B that 
does not change in time (other than with stochastic noise). 

Reweighting. The waiting time to obtain a steady-state 
constrains the efficiency of obtaining a MFPT by measuring 
fluxes via equation [T] This waiting time can vary from the 
relatively short time scale of intra-state equilibration for sim- 
ple systems, to much longer time-scales, on the order of the 
MFPT itself for more complicated systems. To reduce this 
waiting time, we use the steady-state reweighting procedure 
of Bhatt et al. J46). This method measures the fluxes between 
bins to obtain a rate-matrix for transitions between bins, and 
uses a Markov formulation to infer a steady-state distribution 
from the (noisy) data available. 

For instance, let fw,} be the set of bin weights (i.e the sum 
of the weights of the trajectories in each bin), and let [wf] be 
the set of steady-state values of the bin weights. If fy is the 
flux of weight into bin ; from bin j, then in steady-state, since 
the flux out of a bin is equal to the flux into it, 

dw" 

~± = £ i/u -fji) = Z - k fi <) = • (2) 



Since the flux of weight into bin i from bin j is the prod- 
uct of a (constant) rate and the (current) weight in a bin, i.e. 
fij = kjj Wj (true for both steady state and not), we can use 
Eq. [2] to find the inter-bin rates. By measuring the inter-bin 
fluxes and the bin weights, we can approximately infer the 
transition rates, and then find a set of weights that satisfy Eq. 
[2] Once the set of bin weights is found, the weights of the in- 
dividual trajectories in the bins are rescaled commensurately. 

The steady-state distribution of weights thus inferred is not 
necessarily the true steady-state of the system, but it tends to 
be closer to it, and an iterative application of this procedure 
can converge to the true distribution fairly rapidly. In practice 
it has been shown to accelerate the system's evolution to a true 
steady-state by orders of magnitude in some cases [46 1 . 



C. Estimation of Computational Efficiency 

Since it is important to assess new approaches quanti- 
tatively, we compare the speedup in computing time from 
weighted ensemble to a brute-force simulation, (i.e. SSA). For 
a given observable (e.g., the fraction of probability in a spec- 
ified tail of the distribution) and a desired precision, we esti- 
mate the efficiency using the ratio: 



E := 



dynamics time in brute-force SSA 
dynamics time in WE-SSA 



(3) 



Since both WE and brute-force use the same dynamics en- 
gine/software, we can estimate the speedup of WE over brute- 
force by just keeping track of how much total "dynamics time" 
was simulated in each. We employ this measure when estimat- 



ing the advantage of using WE to investigate the tails of prob- 
ability distributions, as well as for finding MFPTs in bistable 
systems. 

Another measure of efficiency we employ for MFPT esti- 
mation gauges how fast WE attains a result that is within 50% 
of the true result (determined from exact or extensive brute- 
force calculation): 



£50% : - 



dynamics time in brute-force SSA to get ± 50% exact 

dynamics time in WE-SSA to get ± 50% exact 

(4) 

This is an assessment of how well WE can extract rough es- 
timates of long time-scale behavior from simulations that are 
much shorter than those timescales. 

Brute-force SSA simulations can be run for long times 
without seeing a transition from one macro-state to another. 
To take account of the brute-force simulations where no tran- 
sitions occurred we use a maximum likelihood estimator for 
the transition time, based on an exponential distribution of 
waiting times, which is a valid approximation for the one- 
dimensional and two-state systems studied below: 



1 n 



Umle 



(5) 



where T is the length of the brute-force simulations, N is the 
number of these simulations performed, n is the number of 
these simulations in which a transition from one state to an- 
other is observed, and f, are the times at which the transition 
is observed. 



D. Limitations of Our Implementation 

We used two different implementations of the weighted 
ensemble framework: WESTPA, written in Python, is the 
most feature-rich and stable If50l . which will be available at 
|http : //chong . chem.pitt . edu/WESTPA| Another, writ- 
ten by Bin Zhang H4l and modified by us, is written in C, 
and is faster though less robust, and is available at |http : | 
|/ /donovanr . github . com/WE_git_code| 

Weighted ensemble (WE), as a scripting-level approach, in- 
herently adds some unavoidable overhead to the runtime of the 
dynamics. This overhead, in theory, is quite minimal: stop- 
ping, starting, merging, and splitting trajectories are not com- 
putationally costly operations. A key issue in practical imple- 
mentations, though, is how long the algorithm actually takes 
to run, i.e the wall-clock running-time for dynamics (here, the 
SSA). 

In practice, overhead can be significant for very simple sys- 
tems, for the sole reason that reading and writing to disk takes 
so much time compared to how long it takes to run the dy- 
namics of small models. In our implementation, data is passed 
from the dynamics engine to WE by reading and writing files 
to disk. This handicap is an artifact of our interface, which 
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could, with minimal work, be modified to something more ef- 
ficient. As a proof-of-principle, the version of WE written in 
C was modified, for the Schlogl reactions and the futile cy- 
cle, to contain hard-coded versions of the Gillespie direct al- 
gorithm for those systems, so as to obviate the I/O between 
WE and BNG. With these modifications, it was difficult to 
ascertain any significant overhead costs at all, and our runs 
completed in a matter of seconds. We also note that as model 
complexity increases and more time is devoted to dynamics, 
the overhead problem becomes negligible. Practical applica- 
tions of WE will, by nature, target models where dynamics 
are expensive, rather than toy models, where they are cheap. 



2. WE Parameters 

The WE data was generated using 101 bins of unit width on 
the coordinate S\. We employed 100 trajectory segments per 
bin that were run for 100 iterations of a t — Is time-step, with 
no reweighting events. The brute-force data is from 10,200 
100-second runs, which is an equivalent amount of dynamics 
to compute as the single WE run, if all the bins were full all 
the time. However, since the bins take some time to fill up, the 
WE run employed only 840,000 one-second segments, which 
makes the comparison to brute-force SSA more than fair. 



III. Models & Results 

We study four different models, ranging in complexity from 
two chemical reactions governing one chemical species, to 
3680 reactions governing 354 species. The models we employ 
are coupled stochastic chemical reactions, which we imple- 
ment and simulate in BioNetGen using the SSA ifTSl [T8l [191 . 
As depicted in Fig.[T| these simulations are, in turn, managed 
by a weighted ensemble procedure. 



A. Enzymatic Futile Cycle 
1 . Model 

The enzymatic futile cycle is a simple and robust model that 
can, in certain parameter regimes, exhibit qualitatively differ- 
ent behavior due to stochastic noise |59, 60]. This signaling 
motif can be seen in biological systems including GTPase cy- 
cles, MAPK cascades, and glucose mobilization |59, 61 , 62]. 

The enzymatic futile cycle studied here is modeled by: 

*/ k, 
Ei+Si ^Bi -+Ei + S 2 

*/ 

(6) 

Kf 

E2 + Si ^ B2 — > £2 + S 1 

*/ 

where kf = 1.0 and k s — 0.1. Here S 1 can bind to its enzyme 
Ei, and in the bound form, B\, (i.e. B\ — E\ ■ Si), it can be 
converted to S 2, and then dissociate (and similarly for 5 2 — » 
Si). The total amount of substrate, Si + 5 2, is conserved, as 
are the amounts of the different enzymes E\ and E2, of which 
is supplied only one of each kind. Following Kuwahara and 
Mura |29l , in the specific system we look at, we set S 1 + S2 = 
100 and E\ +B X =E 1 + B 1 = \. 

Thus constrained, the above system of reactions can be 
solved by a 404-state chemical master equation (CME), to 
obtain an exact probability density for all times when initial- 
ized from an arbitrary starting point. We start the system at 
S\ — S2 =50 and Ei - Et - 1, and are interested in 
the probability distribution of Si after 100 seconds, that is, 
P(S 1 = x ,t= 100). 



3. Results 

Fig. [2] shows that the brute-force SSA is unable to sample 
values of S 1 much outside the range 30 < S 1 < 70, whereas 
the WE method is able to accurately sample the entire dis- 
tribution. Waiting for the brute-force approach to sample the 
tails would take ~ 1/P(tail) ~ 1/10~ 23 ~ 10 23 brute-force 
runs. With a conservative estimate of ~10 4 runs per second, it 
would take ~10 19 seconds, or many times the age of the uni- 
verse, for brute-force SSA to sample the tails at all. WE takes 
2-3 seconds to sample them (note the comparison to exact dis- 
tribution provided by the CME), for an approximate efficiency 
increase E ~ 10 18 . 

For the sake of clarity, error-bars were omitted from Fig. [2] 
Over most of the data range, the error is too small to see on the 
plot. In the tails (of both SSA and WE-SSA) the error is not 
computable from a single run, since there are plot points com- 
prised of only a single trajectory. The error in the estimate of 
the distribution can be inferred visually from the data's depar- 
ture form the CME exact solution. For SSA, however, gen- 
erating uncertainties far all values is essentially impossible. 
When computing quantitative observables reported below, we 
employ multiple independent runs to procure standard errors 
in our estimate. 

From the distribution, we are able to read off useful statis- 
tics. Instead of computing MFPTs, Kuwahara and Mura 11291 . 
and Petzold and Gillespie et al. [30] defined a related quantity, 
the probability of a system to pass from one state to another 
in a certain time: P(x,- — > Xf\At) [29-33]. For instance, one 
might desire to know the probability of the futile cycle to have 
a value of S 1 > 90 at Af = 100. Since WE gives an accurate 
estimate of P(x, t), all that is required to find such statistics 
is to sum up the area under the state of interest. For the fu- 
tile cycle, we find a value of 2.47 x 10~ 18 ± 3.4 x 10~ 19 at 
one standard error for P(S 1 > 90, t = 100 s), as computed 
from ten replicates of the single WE run plotted in Fig. [2] The 
CME equation gives an exact value of 2.72 x 10~ 18 . Sam- 
pling this tail of the distribution at all by brute-force would 
take ~ lAP(tail) ~ 1/(2.72 x 10~ 18 ) ~4x 10 17 brute-force 
runs. Using ten replicate runs, WE is able to sample it us- 
ing 8,400,000 WE segments, which is equivalent to 83,170 
brute-force trajectories, resulting in an increase in sampling 
efficiency by afactor of E ~4xl0 17 /83, 170 ~ xlO 12 for this 
observable. 
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FIG. 2. The probability distribution of 5 1 in the Enzymatic Futile 
Cycle System, after t = 100 seconds, when initialized from a delta 
function at 5 1 = 50, E t = E 2 = 1 at t = 0. The exact solution, pro- 
cured via the chemical master equation (CME), is compared to data 
obtained using the SSA in a weighted ensemble run (WE-SSA), and 
to ordinary SSA, when each are given equal computation time. WE 
data is from a single run. Error bars are not plotted; for a discussion 
of uncertainties, see Sec. |III A3\ 



B. Schlogl Reactions 



1 . Model 



The Schlogl reactions are a classic toy-model for bench- 
marking stochastic simulations of bistable systems [63 1. They 
are two coupled reactions with one dynamic species, X: 



A + 2X ^ 3X 

h 

B^±X 

1-4 



(7) 



where h = 3 x 10"', k 2 = 10 , k 3 = 1Q~\ k 4 = 3.5, A = 
10 5 , and B = 2 x 10 5 . The species A and B are assumed to 
be in abundance, and are held constant. Both the mean first 
passage times and the time-evolution of arbitrary probability 
distributions can be computed exactly ll64l . 



2. WE Parameters 

The WE data in Fig.[3]was generated using 802 bins of unit 
width, 100 trajectory segments per bin, a time-step t = 0.05 s, 
and run for 101 iterations of that time-step, with no reweight- 
ing events. The brute-force data is from 80,200 5-second runs, 
which is an equivalent amount of dynamics as a single WE 
run, if all bins are always full. Were that the case, the WE run 
would compute dynamics for 8,100,200 trajectory segments; 
in our case the WE simulation ran 7,047,300 trajectory seg- 
ments, which makes the comparison to brute-force more than 
fair. 

The WE data in Fig.[4]was generated using 80 bins of width 
10, with 32 trajectory segments per bin, a time-step t = 0.1 s, 



run for 500 iterations of that time-step. Reweighting events 
(see Sec. |II B 2\ were applied every 100, 5, and 2 iterations 
for the data labeled "RW-100", "RW-5", and "RW-2", respec- 
tively. 



3. Results 

Fig. [3] shows how the results of both a brute-force (BF) 
approach, and the WE approach compare to the exact solu- 
tion 11641 . when each employs the same amount of dynamics 
time. We start the Schlogl system with X = 82, i.e. the PDF 
is initially a delta function at X — 82. To investigate rare tran- 
sitions, we study the PDF at time t = 5 s. WE is able to accu- 
rately sample almost the entire distributions, even over the po- 
tential barrier near X = 250, while the BF approach is limited 
to sampling only high probability states. The Schlogl system 
is bistable, with states centered at X = 82 and X - 563, and a 
potential barrier between them, peaked at X — 256. The brute- 
force approach is unable to accurately sample values outside 
of the initial state, and cannot detect bistability in the model. 

For the sake of clarity, error-bars were omitted from Fig. [3] 
Over most of the data range, the error is too small to see on 
the plot. In the tails (of both SSA and WE-SSA) the error is 
not computable from a single run, since there are plot points 
comprised of only a single trajectory. Multiple runs are con- 
sistent with the data shown. The error in the estimate of the 
distribution can be inferred visually from the data's departure 
form the CME exact solution. When computing quantitative 
observables below, we employ multiple independent runs to 
procure standard errors in our estimate. 

WE yields the full, unbiased probability distribution, but 
we again examine an observable investigated by Petzold and 
Gillespie et al. 129143311 . As such, the conversion to the 
rare event statistics of Gillespie et al. is a simple summa- 
tion. From ten replicates of the Schlogl run plotted in Fig. [3] 
the probability that X > 700 at t — 5 seconds, i.e. 
P(X > 700, t = 5 s), is 1.143 x 10~ 9 ± 4.7 x 10" 11 
at 1-cr. The CME exact value is 1.148 x 10~ 9 . Since it would 
take at least 1/(1.15 x 10 9 ) ~ 10 9 brute-force trajectories to 
sample the probability that X > 700 at t — 5, we can estimate 
an improvement in efficiency of using WE over brute-force of 
10 9 /802,000 ~ 10 3 . 

We also estimate the mean first passage time (MFPT) of 
the Schlogl system, which can be computed exactly Il64l . 
Weighted ensemble can estimate the MFPT using Eq.[T|when 
the system is put into a steady-state. For the run that was 
reweighted every 100 iterations, Fig. [4] shows the WE esti- 
mates of the flux from the initial state (X = 82) to the final 
state (X > 563) converge to the exact value in about 100 itera- 
tions of weighted ensemble splittings and mergings, which is 
when the system relaxes from its delta-function initialization 
to a steady-state. The attainment of steady-state is accelerated 



by more frequent reweighting (see Sec. II B 2 on reweighting), 
as is shown in Fig. [4] in the runs that are reweighted every 2 
and 5 iterations. These more frequently reweighted runs yield 
fluxes close to the exact value within about 30 iterations. 
To quantify WE's improvement over brute-force in the es- 
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1 r 




Population of X at t = 5 sees. 

FIG. 3. The probability distribution of X in the Schlogl system, at 
; = 5 seconds, when initialized from a delta function at X = 82. The 
exact solution from the chemical master equation is compared to data 
obtained using the SSA in a weighted ensemble run (WE-SSA), and 
to ordinary SSA. WE data is from a single run. For a discussion of 
uncertainties, see Sec. |IIIB"3l 



timate of the MFPT, we use the measure £50% defined in Sec. 
|IIC| A brute-force estimate of the MFPT would require, op- 
timistically, computing an amount of dynamics on the order 
of the MFPT itself (approximately 5 x 10 4 ) seconds. Since 
transitions in this system follow an exponential distribution, 
the standard deviation of the first passage times is equal to the 
mean of them. WE's estimate of the MFPT is within 50% of 
the exact value after about 30 iterations of WE simulation, at 
which point about 1100 trajectory segments have been prop- 
agated, which is equivalent to propagating about 110 seconds 
of brute-force dynamics. Thus we find £50% ~ 5 x 1 4 / 1 1 ~ 
500. As can be seen in Fig. |4] this value is about a 3-5 fold 
increase over the WE results when reweighting very infre- 
quently (every 100 iterations). 



C. Epigenetic Switch 
1. Model 



This model consists of two genes that repress each other's 
expression. Once expressed, each protein can bind particular 
DNA sites upstream of the gene which codes for the other 
protein, thereby repressing its transcription l65l . If we denote 
the /th protein concentration by g,, the deterministic system is 
described by the equations: 

dgi _ a\ 81 

df 1 + (g 2 /K 2 y t 

(8) 

dg2 _ A2 gl 

df ~ l + igil^r r 

where a\ = 156, a 2 = 30, n = 3, m = 1, K\ = 1, K 2 = 1, 
t — 1 . In our stochastic model, our chemical reactions take 
the form of a birth-death process, the propensity functions of 
which are taken from the above differential equations: 

> gi — > 



where k Q = 1/r, k x {g 2 ) = aj[l + (g 2 /K 2 ) n ], k 2 ( gl ) = a 2 /[l + 
(gilK.Tl 



2. WE Parameters 
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FIG. 4. The flux of probability into the target state (X > 563) for 
the Schlogl system. The exact value is compared to WE results, for 
reweighting periods of every 100, 5, and 2 iterations. The inverse of 
the flux gives the mean first passage time by Eq.[T| 



For this system we implemented 2-dimensional bins: 15 
along g t and 31 along g 2 , for a total of 465 bins. The bins 
along the gi coordinate were of unit width on the interval 
[0, 10], and then of width 10 on the interval [10, 50], with one 
additional bin on [50, 00]. The bins along the g 2 coordinate 
were of unit width on the interval [0,30], with one additional 
bin on [30, 00]. 

The WE data in Fig. [5] was generated using 16 trajectory 
segments per bin, a time-step r = 0. 1 s, and run for 500 itera- 
tions of that time-step, with reweighting events applied every 
100 iterations. Fig.[5]shows six independent simulations using 
these parameters, as well as MLE statistics from our brute- 
force computations. Were all the bins full at all iterations, WE 
would compute, for each of the six runs, 3,720,000 trajectory 
segments of length 0. 1 seconds each, which is equivalent in 
cost to running 372,000 seconds of brute-force dynamics. In 
our case, most of the bins never get populated; we computed 
dynamics for 148,855, 149,516, 148,940, 147,351, 146,804, 
and 149,765 segments in the six different runs. In toto, this is 
equivalent to 89,123.1 seconds of brute-force dynamics. 
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3. Results 

Even the state-space of this two-species stochastic system is 
too large to solve exactly, necessitating the use of brute-force 
simulation as a baseline comparison. A brute-force computa- 
tion was performed using the SSA as implemented in BNG. 
753 simulations of 10 6 seconds each were run, and using an 
exponential distribution of MFPTs, the MLE (see Eq. [5]l of 
the mean and standard error of the mean, /^mle and cry,, were 
found to be 1.3 x 10 6 seconds and 6.5 x 10 4 seconds respec- 
tively. 

The WE results are plotted against the brute-force values in 
Fig. [5] where we have used the relation MFPT = 1/flux (Eq.[T} 
to plot the steady-state flux that brute-force predicts. We plot 
the net flux entering the target state as the simulation pro- 
gresses, because this is what WE measures directly; we can 
infer the MFPT using the above relation. Taking the mean of 
each of the six WE runs after the simulation is in steady-state 
(we discard the first 100 iterations), and treating each of these 
means as an independent data point, WE gives a combined 
estimate for the MFPT of 1.3 x 10 6 + 3 x 10 4 seconds at 1-cr. 

WE is able to find an estimate of the MFPT with greater 
precision than brute-force, using the equivalent of 89,123.1 
seconds of brute-force dynamics. The brute-force estimate 
uses 753 x 10 6 seconds of dynamics, yielding a speedup by a 
factor of E ~ 10 4 when using WE compared to brute-force. 

WE is also able to quickly attain an efficient rough 
estimate of the MFPT. A brute-force estimate of the 
MFPT would require, optimistically, computing an amount 
of dynamics on the order of the MFPT itself (~10 6 s). 
In the six different simulations, WE's estimate of the 
MFPT is within 50% of the brute-force value after 
{52,44,37,40,43,42} iterations of WE simulation, at which 
point {10238,8400,6177,6819,7141,7750} trajectory seg- 
ments have been propagated, which is equivalent to prop- 
agating {1023.8,840.0,617.7,681.9,714.1,775.0} seconds of 
brute-force dynamics, the mean of which is approximately 
775. Thus we find a mean £ 5 o% * 1.3 X 10 6 /775 ~ 1725. 
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FIG. 5. Measurements of probability flux into the target state for 
the epigenetic switch system. Six independent WE simulations are 
plotted, as well as the 3-<x confidence interval for the brute-force 
data, which is from 753 trajectories of 10 6 seconds each. The inverse 
of the flux gives the mean first passage time by Eq.[TJ 



D. FceRI-Mediated Signaling 
1. Model 

To demonstrate the flexibility of the WE approach, we ap- 
plied it to a signaling model that is, to our knowledge, con- 
siderably more complex than any other biochemical system 
to which rare event sampling techniques have been applied. 
The reaction network in this model (see supplementary mate- 
rial [URL will be inserted by AIP] fceriji.bngl) contains 354 
chemical species and 3680 chemical reactions |66|. 

This model describes association, dissociation, and phos- 
phorylation reactions among four components: the receptor 
FceRI, a bivalent ligand that aggregates receptors into dimers, 
and the protein tyrosine kinases Lyn and Syk. The model also 
includes dephosphorylation reactions mediated by a pool of 
protein tyrosine phosphatases. These reactions generate a net- 
work of 354 distinct molecular species. The model predicts 
levels of association and phosphorylation of molecular com- 
plexes as they vary with time, ligand concentration, concen- 
trations of signaling components, and genetic modifications 
of the interacting proteins. 



2. WE Parameters 

The WE data in Fig. [6] was generated using 60 bins of unit 
width, 100 trajectory segments per bin, a time-step t = 0.6 s, 
and run for 100 iterations of that time-step, with no reweight- 
ing events. The brute-force data is from 1484 brute-force runs 
of 60 seconds each, which is equivalent to the dynamics time 
employed in attaining the WE data. No attempt was made to 
optimize sampling times or bin widths in WE. 



3. Results 

Fig. [6] shows the probability distribution of activated recep- 
tors in the FceRI-Mediated Signaling model at time t — 60 s. 
The brute-force SSA approach is unable to sample out to like- 
lihoods much below ~10~ 3 , while WE gets clean statistics for 
likelihood values down to ~10~ 15 , for an estimated improve- 
ment in efficiency E ~ 10 12 . 



IV. Discussion 

We applied the weighted ensemble (WE) [43-50 1 approach 
to systems-biology models of stochastic chemical kinetics 
equations, implemented in BioNetGen fP7ll52l . Increases in 
computational efficiency on the order of 10 20 were attained for 
a simple system of biological relevance (the enzymatic futile 
cycle), and on the order of 10 12 for a large systems-biology 
model (FceRI), with 354 species and 3680 reactions. 

WE is easy to understand and implement, statistically ex- 
act B51 . and easy to parallelize. It can yield long-timescale 
information such as mean first passage times (MFPTs) from 
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FIG. 6. Comparison of WE and SSA for the FceRI signaling model, 
which has 354 reactions and 3680 chemical species. The probability 
distribution is shown for the system reaching a specified level of Syk 
activation (the output of the model, which is a sum of 164 species 
concentrations) within one minute of system time after stimulation. 
Results of 1484 SSA simulations of one minute duration are com- 
pared with WE results generated with an equivalent computational 
effort (several CPU hours in each case). 



simulations of much shorter length. As in prior molecular 
simulations [44, 46], WE has been demonstrated to increase 
computational efficiency by orders of magnitude for models 
of non-trivial complexity, and offers perfect (linear) parallel 
scaling. It appears that WE holds significant promise as a tool 
for the investigation of complex stochastic systems. 

Nevertheless, a number of additional points, including lim- 
itations of WE and related procedures, merit further discus- 
sion. 



A. Strengths of WE 

Beyond the efficiency observed for the systems stud- 
ied here, the WE approach has other significant strengths. 
Weighted ensemble is easy to implement: it examines tra- 
jectories at fixed time-intervals, and its implementation as 
scripting-level code makes it amenable to using any stochas- 
tic dynamics engine to propagate trajectories. WE also par- 
allelizes well, and can take advantage of multiple cores on a 
single machine, or across many machines on a cluster; Zwier 
et al. have successfully performed a WE computation on more 
than 1,000 cores on the Ranger supercomputer 1 50 1 . Addi- 
tionally, WE trajectories are unbiased and follow the natural 
dynamics of the system. WE also yields full probability dis- 
tributions, and can find mean first passage times (MFPTs) and 
equilibrium properties of systems. 



into different regions, and are able to merge and split trajec- 
tories so as to enhance the sampling of rare regions of state- 
space. The approaches described differ slightly in the way the 
splitting and merging of trajectories is performed. WE also 
differs from FFS in that WE does not have to catch trajecto- 
ries in the act of crossing a bin boundary; instead WE checks, 
at a prescribed time step, in which bin a trajectory resides. 
This can be advantageous in that no low-level interaction with 
the dynamics engine/software is required in WE. 

The central hurdle to improving efficiency using acceler- 
ating sampling techniques such as WE, FFS, and NEUS, is 
to adequately divide that state-space by selecting reaction co- 
ordinates that are both important to the dynamics of interest, 
and that are slowly sampled by brute-force approaches. Opti- 
mally and automatically dividing and binning the state-space 
is, to our knowledge, an open problem, and one that, for com- 
plex systems, where a target state is unknown, is not always a 
straightforward one to solve, though adaptive strategies have 
been suggested 1 45 , 48 , 67 1 . 

The wSSA approach [29] differs from the above ap- 
proaches. It does not use a state-space approach, but rather 
uses importance sampling to bias and then unbias the dynam- 
ics. WE seems to have comparable performance to wSSA for 
systems to which both can be applied. Since wSSA biases/un- 
biases reaction rates, while WE divides state-space, the advan- 
tage of one over the other may be situation-dependent. The 
ease of implementation of the WE framework would appear 
to scale better with model complexity than current versions of 
wSSA, though for very small models wSSA may outperform 
WE in measuring select observables. 

A limitation common to accelerated sampling techniques 
used to estimate non-equilibrium observables is the system- 
intrinsic timescale: or the "event duration" time ||68l . 
This timescale represents the time it takes for realistic trajec- 
tories to "walk" from one state to another, excluding the wait- 
ing time prior to the event. The event duration is often only 
a fraction of the MFPT, since it is the likelihood of walking 
this path that is low; the time to actually walk the path is of- 
ten quite moderate. That is, the waiting time in an initially 
metastable state can greatly exceed ?/,. WE excels at over- 
coming the low likelihood of a transition, but no accelerated 
sampling technique can overcome %. 

Finally, it should be noted that all state-space methods that 
branch trajectories, including WE, produce correlated trajec- 
tories, due to the splitting/merging events. While such corre- 
lations do not appear to have impeded the application of WE 
to the systems investigated here, future work will aim to quan- 
tify their effects and reduce their potential impact. The present 
work accounted for correlations by analyzing multiple fully 
independent WE runs. 



B. Comparison to Other Approaches 

WE is most similar in spirit to recent versions of forward 
flux sampling (FFS) [36 1 and non-equilibrium umbrella sam- 
pling (NEUS) |[37l . All of these methods divide up state-space 



C. Future Applications 

Beyond potential applications to more complex stochastic 
chemical kinetics models, the weighted ensemble formalism 
could be applied to spatially heterogeneous systems. WE 
should be able to accelerate the sampling of models such as 
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those generated by MCell 169H7T1 or Smoldyn IT721 . perhaps 
using three-dimensional spatial bins. 

It may be possible to integrate WE with other methods. We 
note that the state-space dividing approaches of a number of 
methods (forward flux E0ll34l - [36l . non-equilibrium umbrella 
sampling 11371 [38l . and weighted ensemble [43—50 1), since 
they are dynamics-agnostic, could be combined with other 
methods that accelerate the dynamics engine itself, such as 
the T-leaping modification of Gillespie's SSA and its many 
variants and improvements ll73Tl76l . to yield multiplicative in- 
creases in runtime speedup. 

More speculatively, WE could be combined with parallel 
tempering methods IT7714791 . WE accelerates the exploration 
of the free-energy landscape at a given temperature, and since 



it does not bias dynamics, the trajectories it propagates could 
be suitable for replica exchange schemes. 

For complex models where exploring the state-space via 
brute-force is prohibitively expensive, WE could also be em- 
ployed to search for bistability, or in a model-checking capac- 
ity [80-82 1 to search for pathological states. 
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